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ABSTRACT 

Non-linear evolution of cosmological energy density fluctuations triggers deviations 
from Gaussianity in the temperature distribution of the cosmic microwave background. 
A method to estimate these deviations is proposed. N-body simulations - in a ACDM 
cosmology - are used to simulate the strongly non-linear evolution of cosmological 
structures. It is proved that these simulations can be combined with the potential 
approximation to calculate the statistical moments of the CMB anisotropies produced 
by non-linear gravity. Some of these moments are computed and the resulting values 
arc different from those corresponding to Gaussianity. 

Key words: cosmic microwave background — cosmology: theory — large-scale struc- 
ture of the universe. 



1 INTRODUCTION 

Many experiments have been designed to analyze the Cos- 
. mic Microwave Background (CMB) anisotropy. The signal 
to be detected is the superposition of various components. 
, In the inflationary models considered here, the components 
generated during linear evolution (Sachs- Wolfe, Doppler, 
and so on) are dominant and Gaussian; nevertheless, there 
are subdominant non Gaussian components as those due 
to non-linear gravity (Rees-Sciama), lens distortions, and 
contaminant foregrounds. The characterization of all these 
subdominant components, including deviations with respect 
to Gaussianity, is necessary for data analysis in experiments 
as MAP and PLANCK. The Rees-Sciama effect was stud- 
ied, using N-body simulations, by Seljak (1996) and Tuluie, 
Laguna & Anninos (1996); nevertheless, deviations with re- 
spect to Gaussianity were not considered at all by these 
authors. The third order moment of the Rees-Sciama effect 
was estimated, using second order perturbation theory (not 
a fully non-linear method describing strongly non-linear evo- 
lution), by Mollerach et al. (1995) and Munshi, Souradeep 
& Starobinski (1995). Here, we are mainly concerned with 
the deviations with respect to Gaussianity produced by fully 
non- linear gravity (Rees-Sciama effect). A new general nu- 
merical method specially designed to estimate these devia- 
tions is proposed. It is based on N-body simulations. Using 
this method, some statistical moments are estimated for the 
first time. 

Quantities Q and Q m are defined as follows: = 



Qb + + f^A and Q m — fib + Qd, where Qb, fid and Qa 
stand for the density parameters corresponding to baryons, 
dark matter, and vacuum, respectively. Quantity h is the 
reduced Hubble constant. We work in the framework of a 
standard inflationary (flat) model with cold dark matter, 
having h = 0.65, Q b = 0.05, fl d = 0.25, and fl A = 0.7. The 
power spectrum has been normalized with erg = 0.93 accord- 
ing with Eke et al.(1996), who estimated as to get cluster 
abundances compatible with observations. In such a model, 
scalar energy density fluctuations are initially Gaussian, and 
they remain Gaussian during linear evolution; nevertheless, 
deviations with respect to Gaussianity arise in the non-linear 
evolution (where Fourier modes mix). 

In this paper, we are concerned with the CMB temper- 
ature distribution. This field can be developed in spherical 
harmonics (AT/T = J^TLo Y^L=-t a e™Yi™)- Taking into 
account the cosmological principle, the CMB temperature is 
assumed to be a realization of a homogeneous and isotropic 
statistical field and, as a consequence of this assumption, 
the m dependence of the bispectrum has -in terms of the 
Wigner-3j symbols- the well known form : 

/ , n ( li ii h \ m 

ii i * 1 J M mi 7772 7773 / 

On account of this formula, we can take mi = 7772 = 7773 = 
to compute Be 1 i 2 e 3 and, then, Eq. (1) allows us to calculate 
{ae imi at 2rn2 ae 3m3 ) for other m values. The m dependence 
of the trispectrum is given explicitly by Hu (2001). 
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If non-linear evolution is not considered at all, namely, 
under the assumption that the CMB anisotropy is only pro- 
duced by linear inflationary fluctuations (which are Gaus- 
sian as a result of the generation mechanism), the ae m co- 
efficients appear to be statistically independent Gaussian 
quantities with zero means and variances (| at m j 2 } = Ce- 
In this situation (see Grishchuk & Martin 1997), all the odd 
moments (mean, bispectrum, and so on) vanish, the spec- 
trum is given by {ae irni ae 2 m 2 } = Ce 1 8e 1 e 2 8 mi ,-m 2 , and the 
higher order even moments can be written as functions of the 
Ce quantities , in particular, the relation (| ai m | 4 > = 3C? 
is satisfied. In order to test that non-Gaussianity develops 
during non-linear evolution, deviations with respect to this 
relation and also with respect to a vanishing bispectrum are 
investigated. 

Along this paper, units are chosen in such a way that 
the speed of light is c = I and, whatever quantity "A" may 
be, A L and Ao stand for the A values on the last scattering 
surface and at present time, respectively. 



2 ANGULAR SCALES, BOXES, AND 
RESOLUTION 

We focus our attention on the anisotropy generated while 
the CMB photons move through the gravitational field pro- 
duced by non-linear cosmological structures. It is worthwhile 
to distinguish two types of these structures: (1) mildly non- 
linear structures having spatial scales larger than those of 
the clusters (superclusters, voids, Great Attractor-like struc- 
tures), and (2) strongly non- linear structures with spatial 
scales either smaller or equal to those of the clusters. In case 
(I), computations can be performed by using second order 
perturbation theory as in Mollerach et al. 1995 and Munshi, 
Souradeep & Starobinski 1995 or, perhaps, the Zel'dovich 
approach and its generalizations; nevertheless, case (2) re- 
quires numerical computations based on N-body schemes. 
We focus our attention on strongly non-linear structures 
as clusters, which could produce significant deviations from 
gaussianity, in particular, during the process of violent re- 
laxation. These structures started to be non-linear (to un- 
dergone violent relaxation) at low redshifts of about z < 30 
(z < 3). 

In our flat ACDM model, the potential approximation 
(Martinez-Gonzalez, Sanz & Silk 1990, 1994; Sanz et al. 
1996), can be used to get the following basic equation: 



AT 
T 

and 



--</> L -n-v L 



■ • dx , 



(2) 



A<j> = AirGSa 2 p B , (3) 

where AZ i s the relative temperature variation -with respect 
to the background temperature- along the direction n, and 
the integral is to be computed, from emission at the last 
scattering surface (L) to observation (O), along the back- 
ground null geodesies. Symbols x l (x), 4>, v, p B , 8, a, t, G, 
stand for the comoving coordinates, the peculiar gravita- 
tional potential, the peculiar velocity, the background mass 
density, the density contrast, the scale factor, the cosmolog- 
ical time, and the gravitational constant, respectively. Since 
we are interested in a gravitational effect, the subdominant 



baryonic component can be neglected and, consequently, the 
density contrast - appearing in Eq. (3)- can be evolved us- 
ing N-body simulations. Our simulations start at z = 50 to 
be sure that all the period of non-linear evolution is taken 
into account. The Rees-Sciama effect produced by non-linear 
structures located at redshift z < 50 is given by the third 
term of the r.h.s. of Eq. (2), which can be approximated as 

follows: [4^1 ~ -2 [ to Vcp-dx ~ 2 f to where t 50 

L i J j J t 50 j t 50 at 

stands for the cosmological time at redshift z = 50. Since 
the anisotropy under consideration can be calculated as an 
integral of d(f>(x,t)/dt, the angular scales of this anisotropy 
are not those subtended by the density contrast itself, but 
those corresponding to the potential <f> (see Fullana & Saez 
2000). 

The remainder of this section is a qualitative discussion 
about angular scales and resolution. No accurate computa- 
tions are necessary and spherically symmetric structures can 
be considered. According to the cosmological principle, any 
spherical cosmological structure is compensated at a certain 
distance or compensation radius r = f. The peculiar mass 
(excess or defect) inside a sphere of radius r, hereafter M(r), 
vanishes at r > £ and, taking into account the relation 

w 

which follows from Eq. (3) in the spherically symmetric case, 
the potential <j> does not depend on r for r > £, furthermore, 
this constant potential must vanishes in order to get a good 
asymptotic behaviour. This means that a given compensated 
structure only produces a significant gravitational potential 
in the region < r < £; namely, inside a sphere with diam- 
eter 2£. In the linear regime, the peculiar velocity produced 
by a structure is v oc V0. Although this relation does not 
hold in the central non- linear region of a cluster, it roughly 
holds in the outer regions, where compensation takes place 
and the density contrast is low. Hence, the peculiar velocity 
vanishes for r > £, at the same place where the gravitational 
potential vanishes by compensation. 

If the clustering of clusters is neglected - a good ap- 
proach in our qualitative analysis- and the present mean 
separation between neighbouring clusters is Mpc, then, 
the mean distance, at redshift z, is £/(l + z)h Mpc and, con- 
sequently, given two similar neighbouring clusters, each of 
them should become roughly compensated at the centre of 
the segment joining the cluster centres (where peculiar ve- 
locity vanishes); therefore, the radius of the region where 
the potential of a given cluster -at redshift z- is not neg- 
ligible appears to be £ = 5/2(1 + z)h. Let us now consider 
Virgo cluster, which is a standard cluster. It is well known 
that the peculiar velocity produced by the Virgo cluster on 
the Local Group ranges between 200 and 300 Km/s; hence, 
taking into account that the distance from the Local Group 
to the cluster centre is ~ 20 Mpc (~ 13/i _1 Mpc), we can 
conclude that the Virgo cluster has not been fully compen- 
sated at a distance of ~ lSh^ 1 Mpc from its centre. This 
strongly suggests that, effective compensation should take 
place between ( = Wli" 1 Mpc (£ = 20) C = 20/i _1 Mpc 
(£ = 40); this second value of ( seems a bit large, although 
an admissible value in the case of Abell clusters, for which, 
compensation should be produced at distances greater than 
those corresponding to Virgo-like structures. 

In a flat universe with cosmological constant, the angle 
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Figure 1. I scale of the non-linear gravitational anisotropy pro- 
duced by non-linear structures as a function of their redshifts z. 
Four values of the present scale £ -defined in the text- are con- 
sidered. 



subtended by the distance £/(l + z)h Mpc placed at redshift 
2 is 



A6 - 



3000/(2) 



where 



/(*) = 



drj 



PWi + ^ + Qa] 1 / 2 



(5) 



(6) 



using these formulae, we have calculated the £ value corre- 
sponding to the angular scale A6 (I — n/A6) as a func- 
tion of 2. The results are displayed in Fig. 1 for £ = 20 
(solid line), £ = 30 (dotted line), £ = 40 (dashed line), and 
£ = 50 (dashed-dotted line). In this Figure we see that, from 
z ~ 0.5 to 2 ~ 5, where clusters undergo violent relaxation 
and the non-linear gravitational anisotropy should be mainly 
produced, the £ values range in the interval (100,800). For 
2 < 0.5, clusters have relaxed and, for 2 > 5, they are at the 
first stages of their non-linear evolution. 

For the sake of briefness and clearness, only the scale 
£ — 400 (with m = 0) has been considered in our calcu- 
lations. Let us now justify this choice taking into account 
Fig. 1. The method proposed in this paper might be sen- 
sitive to both the degree of non-linearity and the speed of 
the <f> evolution. Strong non-linearities and large values of 
d<p/dt (violent relaxation) could lead to numerical problems. 
Non-linearity requires good spatial resolution in the N-body 
simulation, and fast <f> evolution requires appropriate time 
steps in both the N-body simulation, and the numerical in- 
tegrations along the line of sight (see below). Fig. 1 shows 
that non-linear clusters of different sizes located between 
redshifts 0.5 and 5 (some of them undergoing violent re- 
laxation) are contributing to £ = 400, whereas only small 
clusters located at 2 ~ 5 (well before violent relaxation) are 
contributing to £ = 800; therefore, £ = 400 seems to be a 
good choice in order to test our method taking into account 
both non-linear evolution and violent relaxation. Moreover, 



in the model under consideration, the first acoustic peak is 
located at £ = 200 and, around £ = 400, there is a minimum 
of the linear dominant Gaussian effects. 

Cluster evolution is simulated in boxes having -at 
present time- 128 Mpc per edge with a resolution of ~ 1 
Mpc. This choice is a balance between the need of having a 
reasonable number of clusters inside the box (order ten) and 
the need of following the non-linear evolution of the gravi- 
tational potential. The minimum scale we resolve, ~ 1 Mpc, 
is typically of the order of the cluster virial radius and, con- 
sequently, non-linear enough. The fine structure within the 
core cluster is not resolved, but it does not seem to be nec- 
essary because the effect under consideration is produced by 
the gravitational potential, which smoothes all the substruc- 
ture on the scales we are interested in. The effects produced 
by substructures (£ << 20 and strong clustering) should ap- 
pear on angular scales corresponding to £ » 400, for which 
more resolution would be necessary. 

N-body simulations are performed by using a standard 
particle-mesh code (Hockney & Eastwood 1988), which was 
used, tested, and described in more detail by Quilis, Ibanez 
& Saez (1998). 

We are interested in non-Gaussianity and, consequently, 
our goal is not the computation of the angular power spec- 
trum (Ct coefficients) corresponding to the third term of 
Eq. (2). This spectrum was given by Seljak (1996) for the 
Rees-Sciama effect, and by Fullana & Saez (2000) for the 
Integrated Sachs-Wolfe effect. We are looking for the bis- 
pectrum, trispectrum and so on. This fact facilitates our 
implementation of N-body simulations. In the case £ = 400, 
according to Fig. 1, cluster- like structures (£ < 40) are con- 
tributing to the Rees-Sciama effect inside the redshift inter- 
val (0.5,5), whereas structures as superclusters, the Great 
Attractor, and so on (£ > 50), are contributing for 2 > 10, 
just when they were almost linear objects; hence, these 
structures can only produce a very small almost Gaussian In- 
tegrated effect for £ = 400 (see Fullana & Saez 2000) , and no 
important deviation from Gaussianity are expected; hence, 
even if we partially include superclusters and larger struc- 
tures (small box size), we can get very accurate estimations 
of the deviations with respect to Gaussianity (associated to 
strong non-linearity). 

Since we are mainly interested in the evolution of the 
gravitational potential of clusters and substructures and, 
moreover, the chosen box is much larger than the size of 
the region where the potential of standard clusters is signifi- 
cant, the Fourier transform can be performed in our box (as 
it is done in next Section); namely, no important artifacts 
are expected as a result of the periodicity of an universe 
filled by boxes identical to the chosen one. 



3 DEPARTURES FROM GAUSSIANITY 

The basic equations necessary for the numerical computa- 
tion of the bispectrum and other statistical moments are 
first derived using the potential approximation and, then, a 
numerical method for the implementation of these equations 
-based on N-body simulations- is described and applied. 

At present time, the observer is assumed to be located 
at a certain point with comoving spatial coordinates x' p . 
The equations of the null geodesic passing by point x l p arc: 
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x l = x l p + \{a)n l , (7) 

where the function A(a) can be written as follows: 

A(a) = k{o)H q - 1 , (8) 

with 

db 



(9) 



Since we are interested in the deviations from Gaussianity 
produced by non-linear gravity, the third term of the r.h.s. of 
Eq. (2) is calculated. Using Eq. (7), this term can be written 
as follows: 



AT 



(Zp,n) = A(x p ,n) = -2 f (V0 • n)\{a)dt , (10) 

where the dot stands for a time derivative with respect to 
the cosmological time t. 

The following Fourier expansion 



(x,t) = 



d 3 ke- ikS <j> % {t) 



(2tt)3/2 

can be combined with Eqs. (7) and (10) to get: 
A(x p ,n) = 

^L_j\(t)dt J d a fce-*^(t)(fc.n)e- 



(11) 



iX(t)(k-n) 



(12) 



where k is the comoving wavevector. Using: (i) Equation 
(12), (ii) the expansion in spherical harmonics 

oo +e 

A(x p ,n) = ^ a im (x p )Y em (n) , (13) 
e=o m=-l 

(iii) the identity 

(jfe • n) e - iA(t)(K ' fl) = 

oo +e 
^k^i l+1 f t (\k) Y em (n)Y em (k/k), (14) 

1=0 m=-£ 

where j'i(x) — (d/dx)jt(x), (iv) the orthonormality relation 
for the spherical harmonics and (v) the Fourier counter- 
part of Eq. (3), which reads — ~(DS^)/(ak 2 ), where 
D — |(1 — Q,A)HoaQ, the following basic equations are eas- 
ily obtained: 



ae m (x p ) 



(2 



8 -^i e J d 3 ke-^ 3 p k- 2 Y em (k/k)F e (k) (15) 



where 

Fi{k) = / f(k,w)j' e (w)dw 



(16) 



w = k\(a), and f(k,w) denotes the function a~ 1 8^{a) writ- 
ten in terms of the variable w. Eqs. (15) and (16) have been 
obtained using the Fourier transform in all the space, we 
have not introduced either boxes or N-body simulations yet. 
These elements will be only required in order to numerically 
solve the basic equations (15) and (16). 

Evidently, equation (15) can be seen as a Fourier trans- 
form from the k momentum space to the x p physical space, 
and the function to be transformed depends on function 



2 - 




-4 - 



Figure 2. Function G(z) = /(fco,u>(z)) (see text) for an certain 
node ko of the computational box. Each point corresponds to a 
time step of the N-body simulation. 



F((k) given by Eq. (16) and other well known functions. 
In the physical space, each point is the position of a CMB 
observer who expand the temperature contrast in spherical 
harmonics to get the coefficients ae m (x p ). From a mathe- 
matical point of view (to numerically compute ai m (x p )), 
we can introduce an appropriate box in physical space with 
a large enough size (periodic universe) . There is another box 
in momentum space associated to our x p box. In each node 
of the k box, the function 5g(a) must be evaluated a certain 
number of times to perform the time integral involved in Eq. 
(16). As it is shown below, this integral can be performed 
using only the <5g(a) values calculated at each time step by 
the N-body simulation (see below for details). According to 
our arguments in Section 2, boxes of 128 Mpc size contain- 
ing 128 3 observers are appropriate because (a) the resulting 
function <5g(a) has all the necessary information about non- 
linearity and (b) periodicity is not expected to be relevant. 

Quantities at m {x p ) are evaluated on the nodes, x p , of 
the assumed 3D grid as follows: (i) At each time step, ti, of 
the N-body simulations (see above), quantity <5g(a) is cal- 
culated on the nodes of our 128 3 cubic grid (in momentum 
space). At a given time ti, only the quantities f(k, 
and f(k,u>i) are stored, (ii) the integral in Eq. (16) is solved 
for every node and in each time interval as it is de- 

scribed below and, then, results from all the time intervals 
are added to obtain Fe(k). It has been verified that, inside 
any of the intervals (tj_i, U), function fi(k,w) can be very 
well approximated by a straight line and, as a consequence 
of this fact, the integral in Eq. (16) can be analytically calcu- 
lated between ti-i and ti. The linearity of fi(k,w) is due to 
the fact that the time step of the N-body simulation is short 
enough; in order to display this linearity, we present Fig. 2, 
which shows function G(z) — f(ko,w(z)) for a particular 
node ko- Each point is a time step in the N-body simula- 
tion and the linear behaviour between contiguous points is 
evident. The same has been verified for many nodes. From 
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Fig. 2 it follows that function fi(k,w) also has a linear be- 
haviour between more separated points, this fact could be 
used in future to reduce CPU time; nevertheless, in this pa- 
per, linearization is assumed between neighbouring points. 
The equation of the line joining two neighbouring points 

fi(k,w) = wbi(k) + a(k) (17) 

involves two coefficients to be computed from the two stored 
values f(k,Wi-i) and f(k,Wi). Then, using Eq (17), the 
value of the integral (16) can be analytically found in any 
interval [£j_i,tj], the result is 

F e (k) = 

^ [wibi(k) ± a(k)] ji(wi) - [wi-ibi(k) + Ci(k)] je(w t -i) 

i 

-bi(k) / j e (w)dw ; (18) 

w i — l 

it is noticeable that this formula only involves: the coeffi- 
cients appearing in Eq. (17), the spherical bessel functions, 
and the integral of these functions appearing in the last 
term, and (hi) after computing Fi(k), Eq. (15) -a Fourier 
transform- is used to obtain a em (x p ) on the 128 3 nodes of 
the i p grid. Thus, the expansion in spherical harmonics of 
the CMB temperature contrasts corresponding to 128 3 ob- 
servers is obtained and, afterwards, these data are used to 
perform averages. 

Taking into account that a4oo,o is real, we calculate the 
i-averages Ml 00 = (cuoo.o), Cl 00 = (atoo.o), B\ 00 = (aloo.o), 
and /Qoo — ( a too,o) over the 128 3 nodes of the i-th simu- 
lation. Since these averages have been performed on a lim- 
ited sample of rather close observers located in a particular 
region of the universe (nodes in our grid), they appear to 
change from simulation to simulation. In order to calculate 
the statistical moments of the CMB, the simulations must 
be repeated a large enough number of times (to have a big 
enough sample of CMB observers located in independent 
regions). Then, the i-averages corresponding to n different 
simulations (i — 1,2, ...,n) must be averaged again to get 
the quantities M™ 00 , C400, BJ 00 , and K2qq, which must tend 
to the required CMB moments M400, C400, B400, and #400, 
respectively, as n tends to infinity. 



4 DISCUSSION AND RESULTS 

In this paper, a new method has been designed with the es- 
sential aim of calculating deviations from Gaussianity due to 
the Rees-Sciama effect. Our method is different from other 
ones previously used to estimate deviations from Gaussian- 
ity due to weak lensing. According to White & Hu (2001), 
there are too much scales contributing to lensing and sim- 
ulating the full range of scales implied is currently a practi- 
cal impossibility; sophisticated techniques as "special projec- 
tions" (Jain, Seljak & White 2000) or "tiling" (White & Hu 
2001) have been designed to circumvent the problem; never- 
theless, we are not considering lensing here and, as discussed 
in Section 2, the estimate of the deviations with respect 
to Gaussianity produced by the Rees-Sciama effect only re- 
quires simulations involving a moderated range of scales, 
which can be performed (see Section 2). In addition, in the 
Rees-Sciama case, some analytical calculations (Section 3) 



have given formulae to simultaneously calculate the ai m {x p ) 
quantities of many observers, from which, the successive sta- 
tistical moments are calculated by performing appropriate 
averages. 

Condition (17) plays a fundamental role in our numer- 
ical method. It leads to Eq. (18), which allows us to com- 
pute Fi(k) using a few stored data from the N-body simu- 
lation and the spherical Bessel function ji. This computa- 
tion of Fi(k) is much faster than any other procedure based 
on the direct numerical calculation of the integral in Eq. 
(16). In spite of this fact, the quantities involved in Eq. (18) 
are (in this paper) evaluated step by step at each node of 
the computational grid and, consequently, the calculation of 
a,im(Xp) -although feasible in modern computers- has not 
a low computational cost. It is worthwhile to remark that, 
once the quantities ai m (x p ) have been calculated, we can 
obtain, not only one, but many different moments of the 
temperature distribution with very small effort (by comput- 
ing the necessary averages on nodes and realizations). 

With the essential aim of finding non-Gaussian features 
and proving that the proposed method works, we have only 
considered the case £ = 400, m = 0. For very different i 
values computations would be similar, but the range of spa- 
tial scales, the box, and the resolution should be chosen 
again taking into account Fig. 1. The larger the £ value, 
the greater the necessary resolution. More systematic esti- 
mations for other I values -computationally expensive- are 
in progress. According to some comments in Section 1, the 
m dependence can be obtained from the results of the m=0 
case. 

We have performed 90 simulations, in which, the value 
of I Ml 00 I is typically a few times 10 -19 , whereas about 70% 
(25%) of the absolute values of 0400,0 are of the order of 10 -9 
(10~ 10 ). This means that -after ensuring the required nu- 
merical precision- the positive and negative values of 0400,0 
strongly cancel among them to give a | M400 | value as 
small as reported above. On the contrary, the positive and 
negative values of a| o,o (most of them with an absolute 
value between orders 10~ 27 and 10~ 30 ) do not undergo such 
a strong cancellation, and a significant bispectrum appears 
(see below). 

Quantities CJoo, -Blooi and -K400 1 are displayed in Fig. 
3 (in terms of the variable n, from n = 3 to n = 90). 
As the number of simulations n increases, the values of 
these quantities should approach the cosmological limit val- 
ues C400, B400, and ^400, respectively. Of course, the three 
curves in this panel seem to approach constant values as 
n increases. In order to perform a quantitative study of 
this behaviour, we have considered the n-interval [20,90], in 
which all the curves of Fig. 3 seem to approach their limits 
and, in this interval, we have computed the means and vari- 
ances of the above n-quantities. Thus, we have found that, 
within 2a confidence level, the following relations are satis- 
fied: C400 x 10 18 = 1.04 ±0.05, B400 x 10 29 = -4.32 ±0.27, 
#400 x 10 36 = 2.45 ± 0.26. The mean and variance (from 
n = 20 to n = 90) of the quantity rj 00 = 3{C2 00 ) 2 / K2 00 
lead to the following relation F4oo = 1.34 ± 0.09, which is 
valid within 2a confidence level; therefore, for I = 400 and 
m — 0, a non-vanishing component of the bispectrum and 
a non-Gaussian value of the ratio T4oo have been found. It 
is noticeable the fact that quantities C400, -BJooi and K2 00 
really converge as n increases. There are no uncontrolled er- 



6 Aliaga et al. 



1 



ph/0105117 

Jain B., Seljak U., White S., 2000, ApJ, 530, 547 
Martinez-Gonzalez E., Sanz J.L., Silk J., 1990, ApJ, 355, L5 
Martinez-Gonzalez E., Sanz J.L., Silk J., 1994, ApJ, 436, 1 
Mollerach S., Gangui A., Lucchin F., Matarrese S., 1995, 
ApJ, 453, 1 

Munshi D., Souradeep T., Starobinski A. A., 1995, ApJ, 454, 
552 

Quilis V., Ibanez J.M., Saez D., 1998, ApJ, 502, 518 
Sanz J.L., Martinez-Gonzalez, E., Cayon L., Silk J.L., 
Sugiyama N., 1996, ApJ, 467, 485 
Seljak U., 1996, ApJ, 460, 549 

Tuluie R., Laguna P., Anninos P., 1996, ApJ, 463, 15 
White M., Hu W., 2000, ApJ, 537, 1 



5 




36 









" 8 



-I 




29 



5 







20 



40 



60 



80 



n 

Figure 3. Quantities C*J 00 X 10 18 , BJ 00 X 10 29 , and K% QQ X 10 36 
vs. the number of simulations n. 

rors producing strong oscillations and avoiding convergence, 
this verification can be seen as a severe test which has been 
successfully passed by the proposed computational method. 

The Rees-Sciama effect was studied by Seljak (1996) in 
the standard CDM model (A = 0). Four cases correspond- 
ing to different values of the parameters as and Q. m oh 2 were 
considered by this author, who obtained C400 values between 
10~ 17 and 10~ 18 corresponding to a signal, (AT/T) rms , be- 
tween 10 -6 and 10~ 7 . Although we have not computed all 
the Ct quantities, the resulting C400 is comparable to the 
minimum values obtained by Seljak; hence, the non-linear 
gravitational effect should be of the order of 10~ 7 , which is 
below the observational sensitivities of current and planned 
experiments. 

The model of structure formation assumed here is sug- 
gested by far supernovae and recent CMB observations. 
However, other models cannot be ruled out yet. Regard- 
less the structure formation model, our method is ready to 
do estimations. In any case, our results are challenging to 
understand the data in forthcoming CMB experiments 
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